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ABSTRACT 

Protostellar feedback, both radiation and bipolar outflows, dramatically affects the fragmentation 
and mass accretion from star-forming cores. We use ORION, an adaptive mesh refinement (AMR) 
gravito-radiation-hydrodynamics code, to simulate low-mass star formation in a turbulent molecular 
cloud in the presence of protostellar feedback. We present results of the first simulations of a star- 
forming cluster that include both radiative transfer and protostellar outflows. We run four simulations 
to isolate the individual effects of radiation feedback and outflow feedback as well as the combination 
of the two. We find that outflows reduce protostellar masses and accretion rates each by a factor of 
three and therefore reduce protostellar luminosities by an order of magnitude. This means that, while 
radiation feedback suppresses fragmentation, outflows render protostellar radiation largely irrelevant 
for low- mass star formation above a mass scale of 0.05 M . We find initial fragmentation of our cloud 
at half the global Jeans length, around 0.1 pc. With insufficient protostellar radiation to stop it, these 
0.1 pc cores fragment repeatedly, forming typically 10 stars each. The accretion rate in these stars 
scales with mass as predicted from core accretion models that include both thermal and turbulent 
motions; the accretion rate does not appear to be consistent with either competitive accretion or 
accretion from an isothermal sphere. We find that protostellar outflows do not significantly affect 
the overall cloud dynamics, in the absence of magnetic fields, due to their small opening angles and 
poor coupling to the dense gas. The outflows reduce the mass from the cores by 2/3, giving a core 
to star efficiency, e core ~ 1/3. The simulations are also able to reproduce many observation of local 
star- forming regions. Our simulation with radiation and outflows reproduces the observed protostellar 
luminosity function. All of the simulations can reproduce observed core mass functions, though we 
find they are sensitive to telescope resolution. We also reproduce the two-point correlation function 
of these observed cores. Lastly, we reproduce IMF itself, including the low-mass end, when outflows 
are included. 



1. INTRODUCTION 

The origin of the stellar initial mass function (IMF) 
is one of the most fundamental problems of star forma- 
tion. The IMF can be described by single power law for 
stars above 0.5 ( Salpeter|1955[ ), and a broken power 
law ( |Kroupa 2002) lor stars below this mass. Alterna- 
tively, it can be described as a log- normal distribution 
with characteristic mass m c = O.2M that j oins with 
the S alpeter power law for stars above 1.0 M ( jChabrier 
2005| . Any theory of the IMF must explain both the 
functional form and the characteristic mass. A tantaliz- 
ing observational clue to the functiona l form lies in dust 
observations in star-forming regions (|Motte et al.|[ l998 ; 
m " f ~ ' ™~ " 1 ' ' 1 """" Motte 



Testi fc Sargent 1998; Johnstone^amOOO , 2001 
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et al.||2001| IBeuther Sc hilke 2004; Stan ke et al.[ 2006 
Alves et al.|2007| |Knoch et al |2008||Sadavoy et al.|2010p 
These dust maps rind many high density concentrations 
that are consistent with prestellar and protostellar cores. 
When the mass of these cores is calculated, the core 
mass function (CMF) has the same functional form as 
the IMF, but with a higher characteristic mass, ranging 
from 0.2 M to 1 M . If each core is converted to a small 
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number of stars with some efficiency, 0.2 < e core < 1.0, 
the IMF can be recreated. The actual conversion from 
observed core masses to stellar masses m ay not be so 
simple due to cores blending in p rojection (Kainulainen 



|et al.| 



perse 



2009b Michel et al. 



Defor e making stars 



n g m p 
( |Myers| 



sma ll cores that dis- 
|2009 



Padoan & Nord- 



hmdl[20TT1) and cores accreting mass over time (IPadoan 



|& lNordlund|2011[ ). Nevertheless, the CMF likely plays a 
strong role in creating the IMF 
The observed CMF provide s support to core accre- 



tion t heories of star formation (Shu 1977 



McKee & Tan 



2003), which start with a bound core and produce a sin- 
gle stellar system. Simulations of turbulence find the 
functional form of the core mass function (log-normal 
plus power law) is the expected outcome of supersonic 



turbu lence (Padoan & Nordlund 2002 Padoan et al. 



2007). Analytic predictions of a turbulent density field 



with self-gravity can also reproduce this functional form 
flHennebelle & Chabrier||2008l 120091). The characteristic 



core mass is then the Jeans mass at some critical density 
and temperature. However, choosing the correct den- 
sity and temperature is problematic. In purely isother- 
mal turbulence, there is no characteristic Jeans mass. 
As objects collapse and the density increases, the Jeans 
mass decreases. There is no transition where this de- 
crease in Jeans mass will stop without additional physics. 
This means the core masses are either infinitely small or 
functions of the global Jeans mass of the host molecu- 
lar cloud. Observations are consistent with a universal 
IMF, however, even over a range of cloud Jeans masses 
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flKroupa|20Q2l |Chabrier||2QQ3| |Bastian et al ||2Q1Q| ). This 
means trie characteristic core mass must be set by lo- 



by 

cal physics, which isothermal turbulence cannot provide. 
Star-forming regions are approximately isothermal be- 
cause the thermal time scales are much shorter than the 
dynamical time scales, but there are ways to break this 
isothermality. 

One approach is to focus on the coup ling between 
gas and dust in star-for ming environments ( Larson| 2005 



|Elmegreen et al. 2QQ8| ). At low densities, gas- dust cou- 



pling is poor and the gas is theoretically slightly sub- 
isothermal (temperature decreases with increasing den- 
sity). At higher densities, gas and dust are well coupled 
and the gas is theoretically slightly super-isothermal. 
This yields a critical density and temperature at the 
transition that can be converted into a Jeans mass. This 
critical density, p ~ 10 -19 g cm -3 , is lower than the den- 
sities of large star-forming regions like Orion, however, 
and unlikely to explain the characteristic core mass in 
these regions. 

One critical mass is the point when dust becomes 



opaque to its own thermal radiation (Low & Lynden- 
At that density, the gas will heat up and 



Bell||1976|) . 

raise the Jeans mass, creating a minimum Jeans mass of 
fragmentation. A barotropic simplification of this effect 
sets the mass in many simulations 



L£JL 



Bate & Bon- 



ne!Il|2QQ5l IBonnell et al 1F2 006; O ffner et al.||l>0(^| |LSate 
2009a[ IHennebelle et al.||2011|). The density o f this tran- 

13 



sition is ex tremely high (~ 10 13 g/cm 3 ) 



et al. 1998[ ) and the resulting Jeans mass ( 
is much lower than the characteristic core mass 



(Masunaga 
' U.UU4M a 



Lynden-Beli||1976| |Whitworth et al.||2QQ7[ ). In addition 



Low 



the barotropic approximation is inaccurate when com- 
pared to simulations that include dust radiation (|Boss 



eT"aTl[2QQQl [Krumholz et al. 2007 a] IQffner et al.||2 QQ9b; 



Bate|2009b||Price fc Bate|2009||Tomida et al |2Q1QD The 

imp ortance of dust radia tion can be seen in Bate (2009b) 
and Price 



usion of 



Bate (2009), who found that the mc. 
radiation significantly suppresses the formation of brown 
dwarfs despite the near absence of protostellar luminosity 
in the simulations. 

The most powerful break from isothermality comes 
from protostellar radiation. Mass ive protostars are ca- 
pable of heating an entire cloud flKrumholz et al.|[2 007a; 
Cunningham et aT7||2Qll[ |Myers et al.||2011| |Krumholz 



et al.||201ip . Low^mass stars do not have the same long 



range influence, but simulations show they can still dra- 
matically reduce fragmentation in the disk and recover 
1 Mc?) characteristic core mass dOffner et al.||2009b 



jKrumholz et al.||2011| ). Protostellar radiation does not 
create a unique critical density, but it does weaken the 



density dependence of the effective Jeans mass (Bate 



2009bj). 

Given a core mass function, there is still the question 
of CMF to IMF efficiency, e cor e- The primary mecha- 
nism for reducing the core mass is protostellar outflows. 
Stars of all masses show bipolar outflows dur ing their 
formation (Richer et al. 2000| Shepherd 2003) and are 
recreated in MHD sim ulation s with sufficient r esolution 
( jCiardi fc HennebellepOlOl |Tomida et al.||2010[ ). These 
outflows remove mass that would other wise accrete onto 
stars, thereby reducing the final mass (M at zner fc Mc-| 
[Kee|[2Q00l |Arce fc Sargent]|2006l [Wang et al.||2010[ ). An- 



TABLE 1 

Table of simulations 



Name 


Thermal Physics 


Winds? 


B 


Barotropic 


No 


BW 


Barotropic 


Yes 


R 


Radiation 


No 


RW 


Radiation 


Yes 



alytical estimates of mass loss from winds can fully ex- 
plain the range of mass loss expected from observations 
0.2 < e core < 1.0 depending on the de tails of the cores 
and the outflows ( [Matzner fc McKee|2000] ) . The outflows 
travel beyond their stars of origin and deposit energy into 
parsec-scale turbulent motions. Evidence suggests that 
molecular clo ud turbulence appears on the scale of the 
entire cloud (lOssenkopf & Mac Low 2002 Brunt et al.l 



2009), so is most likely driven by sources other than pro- 



tostellar outflows. Nonetheless, the dynamics o n parsec 
scales can be significantly altered by outflows (|Nor man 



Silk|198 0; McKee 1989; Li & Nakamura 2006; tfanerjee 
et al. 1120071 Nakamura fc Li||2007| ISwift fc Welch|l2008 



Carroll et ai.||2010| |Arce et al.||2010| |Wang et al.||2010[ ). 

In this paper, we investigate the fragmentation of ; 



paper, we investigate tne tragmentation of a 
parsec-scale molecular cloud into cores and then into 
stars. This requires refinement to capture the fragmen- 
tation and radiative tr ansfer to fragment at the correct 
mass scale, similar to |Offner et al.| ( |2009b[ ). This also 
requires simulation of protostellar outflows, to ca pture 
the CMF to IMF effi ciency, similar to the work of [Cun- 



ningham et al. (2011) for high-mass star formation. The 
goal of this paper is to explain both the observed CMF 
and the IMF while self-consistently finding e cor e- In or- 
der to do this, we have performed the first simulation of 
a star-forming cluster to include both radiative transfer 
and protostellar outflows. 

We describe our numerical method and simulation 
setup in §2. In §3, we report the results of our sim- 
ulations. We discuss the implications of our results on 
star formation theory and compare to observations in §4. 
We summarize our conclusions in §5. 

2. SIMULATIONS 

We perform four primary simulations with nearly iden- 
tical initial conditions but different physics as controlled 
numerical experiments in order to isolate the impor- 
tance of feedback effects. These simulations all include 
hydrodynamics, gravity and basic sink particle physics 
( Krumholz et al.||2004|), but may also include radiation 
( Krum holz et al. 2007a) and/or sink particle outflows 
( Cunningham et al. 2011). The simulations are shown 
in Table [T] Barotropic simulations are labeled with a 
B, simulations with radiation are labeled with R, and 
simulation with protostellar winds are labeled with W. 
Simulation labels can contain multiple letters. 

2.1. Initial Conditions 

All si mulations have the s ame initial conditions, also 
used in |Offner et al.| ( |2009b[ ). The initial gas tempera- 
ture is T g = 10 A, the box length is L = 0.65 pc and 
the average density is p = 4.46 x 10 -20 g cm -3 , corre- 
sponding to Tin = 1-91 x 10 4 cm -3 . The total box mass 
is 185 Mq. For radiative simulations, the radiation tem- 
perature, T r is initialized to 10 K. The radiation energy 
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erg cm 



density is thus E = aT? = 7.56 x 10" 

To obtain the turbulent initial conditions, we begin our 
simulations without self-gravity and apply velocity per- 
turbations to an initia lly constant den sity field using the 



Mac Low| ( |1999| ). These perturba- 



method described in 

tions correspond to a Gaussian random field with a flat 
power spectrum in the range 1 < k < 2. The application 
of these perturbations continues for three cloud cross- 
ing times and then stops. At this point the turbulence 
follows a Burgers power spectrum, P(k) ex & -2 , charac- 
teristic of supersonic hydrodynamic turbulence. The 3D 
turbulent Mach number is M. = 6.6, which gives a 3D 
rms velocity dispersion, <j v = 1.2 km/s. With this Mach 
number the cloud is approximately virialized: 



5a 2 



C^vir — 



GM/R 



1. 



(1) 



This 

(7 ~ 



is slightly above the 
0.7CR/1 p c) 1 / 2 kms- 1 



linewidth-size relation 



(Solomon et al. 



1987 



Heyer & Brunt 2004), and is equivalent to a 



1.2(R/1 pc ) 1 / 2 km s x , which is w ell within the observed 



2009) 



range (e.g. Falgarone et al 

After driving tor three cloud crossing times, we then 
turn off driving, turn on gravity and follow the subse- 
quent gravitational collapse for approximately one global 
free fall time: 



tff 



3tt 
32Gp 



0.315 Myr, 



(2) 



where p is the mean density of the box. The simulations 
with radiation become prohibitively computationally ex- 
pensive at late times and are stopped at t = 0.83 % 
with a total stellar mass of 30 M for simulation R. The 
barotropic simulations are continued tot = 1.05 t& before 
they are stopped. At this time the total stellar mass in 
simulation B is 50 M compared to the total simulation 
mass of 185 M . There is still gas bound to protostars 
totaling 11 M when the simulations end. Our stellar 
mass estimates may therefore be too low by 20%. 

Given our temperature of 10 K, the Jeans length at p 

^ 1/2 

\j =[^) = 0-20 pc, (3) 



and the Jeans mass is 

4^ (Xj 
3 V 2 



M 



j 



2.7 M . 



(4) 



The turbulent Jeans mass, at density p = M 2 p, is 0.4 
M . 

The calculations have a 256 3 base grid with 4 levels of 
refinement by factors of 2, giving an effective resolution 
of 4096 3 . This resolution corresponds to Ax 4 = 32 AU 
at the finest refinement level. 

2.2. Evolution Equations 

We use the parallel adaptive mesh refinement code 
ORION for our simulations. The numerical method is 



pers (|Krumholz et al. 2007a 


Offner et al.| 


2009b| Cun- 


ningham et al. 2011; Krumho] 


z et al.||2011||Myers et al. 


201ip. ORION solves the equations of compressible gas 



dynamics including self-gravity, radiative transfer, pro- 
tostellar outflows, and radiating star particles, all on an 
adaptive grid. Every cell in the grid has four conserved 
quantities: mass density, p, vector momentum density, 
pv, gas energy density, pe, and radiation energy density, 
E. These conserved quantities can be used to calculate 
derived quantities such as velocity, v, and pressure, P. In 
addition to the gas quantities, we evolve point-mass star 
particles, each with a position x^, mass Mj, momentum 
Pi, angular momentum, and luminosity Li. The sub- 
script i refers to the s tar particle number. T he particle 
method is explained in|Krumholz et al.| (120041) (hereafter 



KKM04), with the addition of radiation (Krumholz et al. 
2007a| ) and outflows ( |Cunningham et al.||2Ull[ ). The lull 



set ol evolution equations tor gas and particles is 



|^+V.(pv)+^[J^ = 0, 



d(pe) 
dt 



V • (pw) = - VP - 



(5) 



$>W(r0 - (6) 



M Wii Vw,iW w (ri)€(0i)'Ti), 



V • [(pe + P) v] = pvV(/> - k r p(4ttB - cE) 



pm H J 

^[eKKMOAWfa) 



(7) 



M Wii W w (ri)m)- 



k B T w K 



dt 



E-V 



/_cA 



M7-ir 

V/: 1 j = K P p ( \~B - cE) + (8) 
A(T,) + ^W(r,), 
(9) 



9 



K pm H 

V 2 (/> = -4^G[p + ^M^(r,)], 
1 



Mi 



-M KKM 04, 



M 7/ 



fu 



-M, 



1 + /. 

Pi = PKKM04, 

0i = acos(ri • ji). 



KKM04j 



(10) 

(ii) 

(12) 
(13) 
(14) 

The quantities entering these equations are defined 
below. Equations ([5| and (|6) are the fluid equations 
for mass and momentum, modified to include particles. 
Equations ^ and ^ are the energy equations for gas 
and radiation respectively. The Poisson equation for the 
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gravitational potential 

part icle evolution is given by equations ([10 ) 



is given by equation ( 91). The 
(TTb and 



(12). We use periodic boundary conditions forall gas 
ana particle quantities. 

For the radiative runs, we adopt Marshak boundary 
conditions for the radiation field. This allows radiation 
to escape from the box as it would from a molecular 
cloud. The equation of state for the gas is given by 



pk B T g 
fim H 



(7- 



l)p(e- y 



(15) 



where \i = 2.33 is the mean molecular weight for molec- 
ular gas of Solar composition and 7 is the ratio of spe- 
cific heats. Since most of the simulation domain is too 
cold to rotationally excite molecular hydrogen, we adopt 
7 = 5/3, representing a monatomic ideal gas. The term 
K,pp(4irB — cE) in equations Q and ([8| represents en- 
ergy exchanged between the raaiation field and the dust 
in our gas, with B = caTg/An representing the Planck 
emission function integrated over all frequencies. The 
opacities np and kr are Planck and Rosseland means 
given by the dust opacities from the iron norma l , com - 



posite aggregates dust model of Semenov et al. (12003) 
We assume that the gas and the dust are thermally cou- 
pled. When the gas temperature exceeds the dust de- 
struction temperature, the energy exchange term goes 
to zero and the gas and radiation unrealistically decou- 
ple. To address cooling from gas above the dust destruc- 
tion temperature, we use t he lin e cooling function A(X^) 
from I Cunningham et aL ( 2006 ) . This removes energy 
fro m the gas and adds t hat en ergy to the radiation field 
(see Cunningham et al. (2011) for further details). The 



radiation flux limiter is given by A = (cothi? — 
where R = |V 'E\j KRpE ( |Levermore fc Pomraning|i 981). 
It should be noted that we have excluded the radiation 
pressure and radiation enthalpy advection terms from 
equations ([6|, |7j) and (j8j) t hat ap pear in the analogous 
equation in jKrumholz etaL (2007a). This approximation 
is justified m the forma tion ol low-mass stars, as shown 
in lOffner et al.j d2009b[ ). 

When radiation is neglected, the energy exchange term 
from equation Q disappears, and we close the system of 
equations with a barotropic equation of state for the gas 
pressure: 



P c lo 



(16) 



where c s q = \fk^T^j 1 prriH is the isothermal sound speed 
at temperature Tq = 10K, 7 = 5/3 and p c is the crit- 
ical density. The critical density determines the transi- 
tion from isothermal to adiabatic regimes and we adopt 
2 x 1 Q~ 13 g cm~ 3 to agree w ith the collapse so 



lution from Masunaga et al. (1998) prior to H2 dissoci- 
ation. Simulations that use the barotropic equation of 
state achieve maximum densities of 5 x 10 -15 g cm -3 , 
with an effective temperature at p m ax of 10.8 K. Most of 
the gas in any given simulations is effectively isothermal. 
The particle quantities Mxkmoa^ Pkkmoa an d 
m Equations ([5] - [Tj) represent the sink parti- 
cle accretion rates of mass, momentum and energy from 
the surrounding gas in the absence of winds as given by 
KKM04. The function W represents a window function 



over the 4-cell accretion zone of the particle and W w rep- 
resents the o utflow window function. O utflows are imple- 
mented as in Cunningham et al. (2011), and summarized 
here. 

Our outflow model is specified by the dimensionless 
parameter f wi which sets the mass flux of outflow as 
a fraction of the accretion onto a star, and v Wj i, the 
wind launch speed. The mass fraction in our simula- 
tions is f w = 0.3. The wind speed is set by the Keple- 
rian speed at the surface of the star, Vk,i — y/GMi/r*^ 
where r*,i is the protostellar radius, but is is capped 
at 60 km s _1 for computational speed. Specifically, 
min(v/ C) i, 60 km s -1 ). The velocity cap has a s im- 



ilar effect to the choice of Cunningham et al. (2011) to 
use v Wy i = Vk,i/S for the most massive stars m the cal- 
culation. Wind gas is injected at temperature T w = 10 4 
K. 

The wind is injected over a window function W Wl 



r 2 if 4Ax < r < 8Ax 
otherwise 



(17) 



which represents a shell just outside of the accretion re- 
gion. The normalization constant C\ is computed nu- 
merically to avoid numerical aliasing effects that occur 
from injecting a spherical wind into a Cartesian grid. 

The exact angular distribution of the wind is described 
in the function £. The functional form is taken from 
|Matzner fc McKee| Jl999| > , 



In 



(sin 2 + 0§) 



(18) 



where #0 is a flattening parameter that sets the opening 
angle of the wind. We use the fiducial value of 6q = 0.01. 
Equation (18) is averaged over the solid angle subtended 
by a grid cell at the outer radius of W w . This averaging 
is particularly important near = 0. In addition, £ is 
set to zero for 6 ~ n/2, so that winds are not injected 
directly into the plane of any equatorial disks. 

We update the luminosity of each star, Z^, using the 
protost ellar evolution model described in Off ner et al.| 
(2009b). In this model, 75% of the accretion energy is ra- 
diated away while 25% is nominally used to power winds. 
The energy of winds in our simulations is already deter- 
mined by f w and v Wj i, which is typically 15% of the accre- 
tion energy. The remaining 10% of the accretion energy 
is effectively lost. When winds are not present, we still 
only radiate 75% of the accretion energy for consistency 
across simulations. 

The evolution equations ca n be described as th e fluid 
and radiation equations from|Offner et al.|Q2009b) com- 



bined with the particl e equations and line cooling of |Cun- 
ningham et al. (|2011|), but there is one important mod- 



ification. In the KKM04 sink particle methodology, all 
particles with overlapping accretion zones are merged to- 
gether. This gives an effective merger radius of 8 cells, or 
256 AU at a grid resolution of 32 AU. To limit this effect, 
we changed the merger radius to 4 cells, representing the 
point when a particle is in the accretion zone of another 
particle. This gives an effective merger radius of 128 AU 
at our resolution. Even with this improvement, our par- 
ticle algorithm will unrealistically merge stars that pass 
within 128 AU. To address this, we have implemented a 
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mass limit of m mer ge = 0-05 M , above which stars do 
not merge. This limit is chosen to correspond to the mass 
of second c ollapse in the formation of a star with final 



mass 1 Mq ( |Masunaga et al.|1998| |Masunaga fc Inutsuka 
[2000] ). Second collapse occurs when the protostar's core 
temperature becomes high enough to dissociate molec- 
ular hydrogen. Before second collapse, protostars are 
extended balls of gas with radii of a few AU and have a 
much higher collisional cross-section than main sequence 
stars. After second collapse, stellar mergers should be ex- 
tremely rare and we do not allo w them. This approach is 
also used in |Myers et al.| ( |2Q11[ ) . The mass of second col- 
lapse depends on the accretion history of the protostar 
and is necessarily less than 0.05 for brown dwarfs 
flStamatellos fc Whitworth|2009|[Bate|2011[) . The effects 
of numerical merger suppression is explored in §4.9. 

ORION utilizes a second order Godunov scheme to 
solve the equations of compr essible gas dynamics (P3> 



elove et al. 



[r¥l 

w 



T998j [Klein 1999). These are equations 
(JT]), excluding terms from stars and radiation. The Pois 
son equ ation (|9| is solved using a multi - ffrid iteration 



scheme ( |Truelove et~aLl[l998j |Klein||1999) |Fisher| [2002 ) . 

The flux-limited diffusion radiation equation (|8| and the 
radiation terms in equation Q are solved using t he con- 
servative update scheme from |Krumholz et al. (2007b) 



modified to include t he pse udo-transient continuation of 
jShestakov fc Offner| ( [2008] ). 

We use the Truelcve criterion ( |Truelove et al.||1997[ ) to 
determine the addition of new AMR grids so that the gas 
density in the calculation always satisfies 



P< 



J 2 7TC 2 S 

G(Axi) 2 ' 



(19) 



where Axi is the cell size on level I. We adopt a Jeans 
number of 0.125. In the simulations with radiative trans- 
fer, it is necessary to resolve the spatial gradients in the 
radiation field. Areas of high radiation gradients are near 
accreting stars, which tend to already be refined under 
the Truelove criterion. This is not always true for more 
evolved stars, which have higher luminosities and have 
accreted the dense gas that would trigger refinement. We 
find that the radiation gradients are adequately resolved 
by refining whenever \VE\Vxi/E > 0.25. 

3. RESULTS 

3.1. Large Scale Evolution 

The evolution of the barotropic simulations with and 
without winds is depicted in Figure [I] Figures [2] and [3] 
depict the evolution of the radiative simulations without 
and with winds, respectively. In all simulations, for 
t < 0.4£ff , there are cloud-scale filaments that slowly con- 
tract, allowing 3 turbulent cores of width ~ 20, 000 AU 
and density ~ 10 -19 g cm -3 to form. This length is half 
the Jeans length at the average cloud density. The first 
core to form is at the center of each panel in Figures [TJ3J 
the second core is left of the central core, near the left 
edge, and the third core is near the right edge, at the end 
of a simulation-wide filament. At this point, the cores be- 
gin to fragment, while new cores form, eventually forming 
6 fully developed cores. These cores each have a cen- 
tral stellar system. In simulation R, this central system 
is all that forms in each core. In all other simulations, 
the central system represents ~ 75% of the stellar mass 




le-01 



le-01 



Fig. 1. — Column density of the entire simulation domain for BW 
(left) and B (right) at times 0, 0.2, 0.4, 0.6, 0.8 and 1.0 % from top 
to bottom. Star particles are marked with white circles. There is 
very little difference on the domain scale with and without winds 
for the barotropic simulations. The color bar is g cm -2 and the 
entire domain is 0.67 pc across. 



in the core and an additional group of low-mass stars 
forms, totaling ~ 10 stars per core. These cores with 
multiple stars generally resemble observed high-stellar - 
density cores ( |Teixeira et al.||2007| |Chen fc Arce||2010[) . 
There are an additional 20 stars m cores that are still 
forming at the end of the simulation, giving a global to- 
tal of 80 stars. Three of the cores coalesce by the end of 
the simulation to form a single group of 30 stars. 

The evolution of the 3D rms velocity dispersion, o~ v , is 
shown in Figure [4] The global turbulence decays un- 
til star formation ramps up at t ~ OMf?. There are 
two main mechanisms for star formation to increase cr v . 
First, as stars accrete mass and deepen their gravita- 
tional potential, the surrounding gas can convert gravita- 
tional energy into kinetic energy as it falls into the stars. 
This is shown in the gradual increase in Mach number 
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Fig. 2. — Column density (left) and density weighted temperature 
(right) for simulation R at times 0, 0.2, 0.4, 0.6, and 0.8 % from 
top to bottom. Star particles are marked with white circles. 

for t > 0.6£ff in the B simulation. This effect is strong 
enough to return a v to near its original virialized value 
by itself. In rare cases, a many-body close encounter be- 
tween stars will eject some gas at high velocities. There 
is not much momentum injected this way and the en- 
ergy quickly dissipates, but it causes spikes in a v for the 
barotropic simulations, which have more small-scale frag- 
mentation and therefore more many-body close encoun- 
ters. The second mechanism occurs when protostellar 
winds are included. Some mass accreted onto stars is di- 
rectly injected around the stars at high velocities. This 
causes the smooth increase in <j v for simulations BW and 
RW as well as spikes from events with particularly high 
accretion rates that lead to bursts in wind momentum. 

The total momentum injected by winds for model BW 
is shown in Figure |5| For comparison, a characteris- 
tic value of the magnitude of the momentum associated 
with internal motions in the cloud, M c \ on &(j v is also plot- 
ted. For t > 0.8£ff, the total momentum that has been 
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Fig. 3. — Same as Figure [2] but for simulation RW. The high 
temperature regions are the paths of outflows. It only takes a 
small amount of gas at 10 4 K to move the average temperature 
above 30 K through that line of sight. 

injected by the winds is greater than the characteris- 
tic cloud momentum. At this point, the total amount 
of turbulent momentum that has been dissipated (in- 
cluding dissipation of wind momentum) is roughly the 
total amount of momentum that has been injected by 
the winds. By the end of the BW and RW simulations, 
the total wind momentum injected into the cloud is over 
twice the characteristic cloud momentum. The kinetic 
energy injected from the winds dissipates over time, sug- 
gesting a steady-state solution where the velocity dis- 
persion of the cloud is constant with time as the winds 
replenish energy as quickly as it can dissipate. 

3.2. Evolution of the Protostellar Population 

The total mass in stars and the total number of stars 
as functions of time is shown are Figures [6]and[7| The 
realization of the initial turbulence is slightly different be- 
tween the radiative and barotropic simulations, so star 
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Fig. 4. — Time evolution of global rms Mach number for simula- 
tions with and without winds and with and without radiation. The 
turbulent energy in all simulations decays for half a global free fall 
time, at which point gravitational potential energy from stars is 
converted into kinetic energy, which raises the rms velocity. When 
winds are included, they contribute over twice as much energy as 
gravity itself. 



however, because protostellar luminosity inhibits frag- 
mentation and the winds reduce that luminosity. The 
three simulations other than R show a dramatic increase 
in the number of stars at 0.6 t& < t < 0.8 




Fig. 6. — Total mass in stars as a function of time for the four 
simulations. The mass of the entire simulation domain is 180 Mq. 
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Fig. 5. — Total momentum that has been injected by winds over 
time for the barotropic simulations. For comparison, the total mass 
of the simulation multiplied by the velocity dispersion is also plot- 
ted. The total wind momentum integrates all injected momentum 
over time, even from winds that have decayed. As a result, the 
amount of injected momentum is eventually higher than the actual 
momentum of the cloud. 

formation begins at different times. This is the result 
of a slightly higher maximum density due to turbulence 
and it is unrelated to the choice of radiative or barotropic 
thermodynamics. The first stars form at t = 0.2% for the 
radiative simulations and t = 0.35% for the barotropic 
simulations. The turbulence overlaps enough between 
the radiative and barotropic cases, however, that at later 
times the total mass in stars at any given time is similar 
for the two cases when winds are not included. Winds 
reduce the mass in stars by about a factor of 3 in both 
the radiative and barotropic cases. The number of stars 
does not change between BW and B, implying winds do 
not cause or suppress fragmentation by themselves. The 
number of stars in RW is significantly greater than in R, 
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Fig. 7. — Number of stars as a function of time for the four sim- 
ulations. In the barotropic case, the number of stars is unaffected 
by winds. In the radiative case, radiation suppresses the number 
of stars unless winds are present. 
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The evolution of median stellar mass is shown in Fig- 
ure [8j This is a rough proxy for the characteristic mass of 
the protostellar mass function. Note that it will always 
be lower than the median mass of the IMF, because not 
all of the stars have finished accreting. Both BW and 
RW maintain a median around 0.05 (similar to that 
of th e protostellar mass functions in Mc Kee fc O ffner 
2010) throughout the simulation. The median does in- 



crease for simulation BW around t > % as the formation 
rate of new stars decreases. The median of B fluctuates 
more, but is around 0.2 M . Lastly, R maintains a me- 
dian around O.5M . This general behavior should be 
expected. The median mass is lowest when winds are in- 
cluded and highest when radiation is allowed to suppress 
fragmentation. The case with both winds and radiation 
ends up similar to BW because winds reduce protostellar 
luminosities. 




Fig. 8. — Median mass of stars as a function of time for the 
four simulations. The two cases with winds maintain low medians 
throughout the simulations. The case with radiation without winds 
(case R) is able to suppress fragmentation and new star formation 
largely stops as the original stars accrete mass. 



The global luminosity evolution for the radiative simu- 
lations is plotted in Figure [9j The winds reduce the total 
luminosity by a factor of up to 10 at any given time. This 
is expected since the global accretion luminosity is 



He 



GM^M+i GM*, tot M* 



-,tot 



(R*) 



(20) 



given that the total mass in stars M* )tot is reduced by 
a factor of 3 when winds are included, and the total 
accretion rate of stars M* 5 tot is also reduced by 3, the 
total luminosity is therefore reduced by a factor of 9 as- 
suming the characteristic stellar radius, (i?*), does not 
change. Main sequence stars typically have a positive 
correlation between mass and radius, suggesting the fac- 
tor of 9 is an upper limit. However, at any given point 
in time, many stars in our simulations are in the de- 
generate regi me where mass an d radius are negatively 
correlated ( |Chabrier et al.|[2009| , which counteracts the 
positive correlation from the higher mass stars and keeps 
the total luminosity ratio near 9. 

The average stellar luminosity is shown in Figure [To} 




Fig. 9. — Total stellar luminosity versus time for simulations with 
radiation, both with and without winds. Winds dramatically lower 
the luminosity. 

As was seen in the plot of total luminosity, the mean 
and median values of protostellar luminosity are much 
lower when winds are included. The disparity in aver- 
age luminosity is even greater than the disparity in to- 
tal luminosity because there are fewer stars when winds 
are excluded. The average luminosity in simulation R is 
heavily influenced by a single 6.6 M star that accounts 
for over half the total luminosity in the simulation. Un- 
like the low-mass stars, most of this luminosity is pow- 
ered by nuclear fusion rather than accretion. Pr otost ellar 
luminosities will be discussed further in section [4731 



„ 120 
c 

£ 60 
ii 40 
h 20 



— R Mean 

- - R Median 










Fig. 10. — Mean and median stellar luminosity versus time for 
simulations with radiation. The top panel is the simulation without 
winds and the bottom panel is the simulation with winds. 



3.3. Thermal Evolution 

All simulations start at a background temperature of 
10 K and are bathed in 10 K radiation. Stellar radiation 
and mechanical energy from protostellar outflows, can 
raise this temperature. We have identified gas heated 
above 12 K as thermally affected by stellar feedback. 
Turbulent dissipation alone heats almost no gas above 
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12 K, leaving stellar feedback as the only explanation for 
this heating. The total mass of this gas is shown in Fig- 
The simulation with winds has significantly less 



ure 
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Fig. 11. — Total gas mass heated above 12 K versus time, com- 
pared to the background value of 10 K. The total mass in stars is 
also plotted for reference. 

heated gas than the simulation without winds. This is 
due to the reduced luminosity caused by the winds shown 
in Figure [9] In each simulation, the mass in heated gas 
roughly follows the mass in stars. When winds are in- 
cluded, the mass in stars drops by a factor of 3, and the 
mass of heated gas also falls due to the reduced stellar 
luminosity. Some of the lost radiative heating is replaced 
by mechanical heating from outflows. This is evident in 
the ratio of heated gas mass to stellar mass. The mean 
value of this ratio is 0.9 without winds, but rises to 1.5 
with them. Wind gas is injected at a temperature of 10 4 
K, but cools quickly. The mass in gas with temperatures 
above 1000 K is less than 2% of the mass in stars at any 
given time. 

To further explore the heated gas, temperature-density 
phase plots with and without winds are shown in Figure 
[12] The phase plots with and without winds are notably 
different in two areas. First, the wind gas fills the high- 
temperature, low-density domain, while the same domain 
is empty without winds. Second, high-density gas with 
p > 10 -16 g/cm 3 has a higher temperature range without 
winds than with winds because the extra stellar luminos- 
ity heats that gas. When winds are included, that dense 
gas is less common in addition to being colder; there is 
more fragmentation, which turns dense gas into stars. In 
addition, some of the gas is also blown away by the winds 
themselves. 

4. DISCUSSION 

4.1. Supporting a Cloud with Outflows 

The turbulent evolution shown in Figures [4] and [5] 
roughly agrees with previous simulati ons of molecular 
clouds wit h outflows, such as those in |Nakamura fc Li| 



(2007k and Wang et al. (2010). Turbulent energy decays 
initially, only to be replaced by kinetic energy from winds 
and from gravity. While these sources can increase the 
total kinetic energy of gas, the new turbulence is funda- 





FlG. 12. — Phase plot showing total gas mass as a function of 
temperature (y-axis) and density (x-axis) for radiative simulations 
with (left) and without (right) winds. Phase plots are taken at 
times of 0.25, 0.5 and 0.75 % from top to bottom. The high- 
temperature, low-density gas on the left part of the wind phase 
plots is outflow gas. Warm, high-density gas is gas near a luminous 
star. 



mentally different from the isotropic, homogeneous hy- 
dro dynaimc_turbujencei^ This result is also 
seen in |Nakamura fc Li| ( |2007[ ), who find that the late 
time turbulent statistics do not match expected isotropic 
hydrodynamic results. One key difference is the energy 
from outflows is highly anisotropic. Outflow cavities are 
marked by long channels with high velocity shear be- 
tween the fast outflow gas and the slow ambient gas. This 
shear is detectable as solenoidal energy. There is some 
compressive energy at the head of the outflow cavity, but 
most of the surface area is the side walls of the cavity 
and not the head. To measure the relative importance 
of solenoidal and compressive motions, we use the ratio 



/V x v\ 2 )/((V • v) 2 ), shown in Figure |13| In the case of 
isotropic turbulence, this ratio is 2, w hicn is also the ratio 
of sol enoidal to compressive energies ( Elmegreen & Scalo 
2004 ). Wind injection greatly increases the solenoidal 
velocities, steadily increasing the solenoidal to compres- 
sive ratio over the course of the simulation. Bursts in 
outflows injection around t = OMf? and t = 0.65£ff are 
also visible in the solenoidal velocity. Anisotropic tur- 
bulence behaves differently than isotropic turbulence; in 
particular, it takes longer to decay since it decays on 
the crossin g time calculated fro m the smallest velocity 



dispersion (Hansen et al. 



20111). It is difficult to mea- 



sure this increase in the decay time in our simulations, 
however, because the winds drive turbulence over a wide 
range of scales. 

The other major difference between the initial turbu- 
lence and wind-driven turbulence is seen in the rms ve- 
locity, Odense? of gas with p > p. Even in isotropic, 
homogenous, hydrodynamic turbulence, there is a neg- 
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Fig. 13. — Ratio of solenoidal to compressive velocities squared 
((|V x v\ 2 )/((V • v) 2 )) versus time with and without winds for 
the radiative simulations. For pure hydrodynamic, isotropic tur- 
bulence, this ratio is around 2. This ratio stays near 2 when winds 
are excluded. When winds are included, the turbulence is much 
more anisotropic, leading to higher solenoidal fractions. 

ative correlat ion between densit y and velocity, causing 
tfdense < o~ v ( Offner et al.||2009a ). The winds themselves 



are collimated, very low density gas and have difficulty 
transmitting energy into high density gas. This means 
that while a v is much greater with winds than without, 
tfdense does not change much when winds are included. 
The evolution of adense compared to a v is shown in Fig- 
ure run 
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Fig. 14. — Time evolution of rms Mach number of dense gas and 
all gas with and without winds. Winds significantly raise the Mach 
number of the light gas, but do not strongly influence the dense 
gas turbulence. 

Because the dense gas is relatively unaffected by the 
outflows, i f our cloud had been ce ntrally concentrated 
like that of |Nakamura fc Li| ( |2QQ7j ) or |Wang et aT] ( |2Q10|| , 
the dense part ol the cloud would have most likely col- 
lapsed on itself even with the support of protostellar out- 
flows. Magnetic f ields may have an effect, as shown by 
| Wang et al.| ( |2010[ ). The magnetic fields help transmit 
outflow energy to a much larger solid angle, so that even 



a centrally concentrated cloud can achieve a quasi-static 
balance between outflows and gravity. 

4.2. Comparison to the IMF 

Our simulations are marginally able to capture bina- 
ries, so most star particles should represent a single stel- 
lar object instead of a stellar system. The typic al binary 
separ ation for main sequence stars is 50 AU ( Ma thieu] 
1994), just slightly larger than our resolution ol 32 AU. 
We cannot form stars within 128 AU of another star 
due to merging star particles, but stars that form be- 
yond 128 AU can spiral in to 32 AU through interaction 
with gas. At any given time, about 2/3 of our star par- 
ticles are in stellar multiples. The multiple properties 
are dynamic due to unstable high-order multiples and 
we cannot compare to the observed system properties. 
We can, however, compare to observed stellar proper- 
ties. The mass functions of the four main simulations 



ar e shown in Figure 15 and compared to the stellar IMF 
in [Chabrier (2005 ) as w e ll as p rotostellar mass functions 
from McKee & Offner (2010b. The mass functions in 
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Fig. 15. — The mass function of all stars in ea ch simulation are 
shown in blue histograms. The stellar IMF from iChabrierl 120051 
is plotted as the solid green line. The protostellar mass functions 
for the turbulent core m odel and the isothermal sphere model from 
|McKee & Offner| pOlO] are the dashed lines and dash dotted lines, 
respectively. Top lert: RW at t = 0.83%. Top right: R at t = 
0.83%. Bottom left: BW at t = 1.09%. Bottom right: B at 
t = 1.03% 



Figure 15 are shown at the latest time available for each 
simulation. The barotropic simulations could be evolved 
to later times due to the computational expense of flux- 
limited diffusion with many stellar sources. The mass 
functions in Figure 15 are not exactly comparable to an 
IMF because some oi the stars are still accreting. The 
barotropic mass functions have evolved to a later time 
and should be similar to the IMF. The radiative simu- 
lations are still actively accret ing and are closer to th e 
protostellar mass functions in |McKee fc Offner (2010). 
Among the protostellar mass functions, the turbulent 



core and competitive accretion (not shown in Figure 15) 
models both roughly match RW. The isothermal sphere 
protostellar mass function is the best fit to mass func- 
tions of the simulations without winds, but this is merely 
a function of both mass functions being too top-heavy. 
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The actual accretion is not described by an isothermal 
sphere, and is discussed further in section [43} 

As expected, the simulations without winds have mass 
functions skewed to higher masses than those with winds. 
The mass functions of both simulations without winds 
have too much mass at the high end compared to the 
IMF. The best fit to the IMF is from the BW simulation. 
It should be noted that the normalized mass functions 
for BW and RW look nearly identical when compared 
at the same time (explored in the next section). It is a 
good assumption that RW would also eventually match 
Chabrier at t ~ % . 

4.3. Comparison to Protostellar Luminosities 

Theoretical predictions of protostellar luminosities are 
often too high compared to observations of regions of 



low-mass star formation feenyon et al. 

J 2UU9[ ), so it is important to 



Evans 2005 Enoch et al 



1990 Young & 



compare our own simulations with observations. The 
mean and median luminosities of protostars observed in 
nearby clusters are (L bs) = 5.3^'g L and L bs,med = 
respectively (Enoch et al. 112009 



+0.7 



0.4 



The 



Evans et al. 
mean and 



1.5 

[2Q09| lOfTner fc McKeej |2011| ). ThcTtypicaT 
median luminosity in our simulation with winds are 
(L) =6.9 Lq and L m ed = 1.4 L , in agreement with the 
observations. An additional useful quantity is the stan- 
dard deviation of the luminosity, cr(L). This is sensitive 
to outliers in the luminosity distribution, but this can 
be mitigated by using the log of the luminosity. We find 
a(logL) = 0.77 dex. This matches the observed value 
cr(logL b s ) = 0.7+01. For these comparisons, we have 
discarded stars with L < 0.05 L , below the detection 
limit of the observations. 

Theoretical frameworks that do match th e observed 
protost ellar lumi nosities are pres ented in |McKee "fc 
|Offnerl ( [20To| and |Offner fc McKee| ( |20TT| ) . These models 
differ from t he st raightfo rward , isothermal sphere mod- 



els (Fletcher & Stahler 1994a b) in three important ways. 
First, lUffner &; McKee| ( |20llp assumed that 1/4 of the 



energy ol the gas that accreted onto the star was re- 
moved non-radiatively. This effect is captured in our 
simulations. Second, they considered accretion rates that 
rise over time, such as predicted in core accretion and 
competitive accretion. These accretion rates all take the 
form, 



rn 



mi 



(21) 



where m is the instantaneous mass of a protostar, rrif 
is the final mass of a protostar, and rai is a constant 
throughout a cloud. More realistic accretion rates will 
rise at early time s and s lowly decline over time, in what 
|McKee fc Offner| ( |2010D call 'tapered accretion 7 There 
are multiple app roaches to this; the one taken by McKee 
& Offner (2010) is to assume a linear decrease in the 
accretion rate with time, m oc 1 — t/tf, which implies 
that 



rn 



m V 



/ m 



i-j" 



1 1/2 



(22) 



The values of j and jf for several different accretion the- 
ories are shown in table [2] as well as the values measured 



in our simulations with winds. Radiative and barotropic 
simulations produce the same j and jf. The fits shown 
are for equation (22). Our data do not agree with the 
functional form for untapered accretion (equation ([21])), 
but if this were used, jf would remain the same while 
j ~ —0.1. Our measured values of j — 0.3 ± 0.2 and 



TABLE 2 

Accretion rate dependencies on instantaneous and final 
protostellar mass 



Accretion Mechanism" 


3 


3f 


Isothermal Sphere 
Turbulent Core Accretion 
Competitive Accretion 
Simulations with Winds 




0.5 
0.67 
0.3 ±0.2 



0.75 
1 

0.6 ±0.2 



a All accretion is tapered, as in equation |22| 

jf = 0.6 ± 0.2 marginally agree with the turbulent core 
accretion model, but do not match any of the other the- 
ories. Our j and jf actually lie in between isothermal 
sphere and turbulent core accretion. This suggests the 
most appropriate theories may be the the two component 
turbulent core (2CTC) model, which is approximately 
isothermal sphere accretion at early times and turbu- 
lent core accretion at later times. Other candidates in- 



clude the TNT model ( [Myers fc Fuller||1992| ), which in- 
cludes both thermal and nonth ermal support, and the 
core-clump accretion model from Myers (2011b) which is 



approximately isothermal sphere accretion at early times 
and reduced Bondi accretion at late times. Our low j 
does not agree with competitive accretion, which accel- 
erates with mass. 

While we a gree with tapered accretion from McKee 



Offner (2010), we do not rule out other tapered models. 



such a s the exponential tapering in Schm eja fc Klessen| 
(2004 |. Our avera g e ac cretion rates are lower than 



|Schmeja fc Klessen (2004) by an order of magnitude 
partially due to our inclusion of protostellar outflows, 
but we see a similar rise and fall of accretion, which can 
be fit by many functional forms. 



T he last effect considered by McKee fc Offner (2010) 



and |Offner & McKee| fl20lT1 ) is episodic accretion from 



FU On type protostars (e.g. [Hartmann fc Kenyon|1996l ) 
If protostars spend most of their lite in a low- accretion, 
low-luminosity phase and accrete nearly all their mass 
during short intense accretion bursts, the median stellar 
luminosity can be greatly reduced. This episo dic accre- 
tion is though t to arise from disk i nstab ilities ( |Vorobyovj 
|fc Basu|[2006l ). |Offner fc McKee| ( |201l[ ) estimated that 
^5% of trie" mass was accreted in this manner. We do not 
resolve disks and therefore do not see episodic accretion 
in our simulations. Given that our simulations match the 
observed luminosity function without episodic accr etion, 
models that rely on substantial episodic accretion ( |Dun-| 
ham et al. 112010 Stamatellos et al.||2011|) may no longer 



be necessary. 

4.4. Suppressing Fragmentation with Radiative 
Feedback 

Fragmentation in the simulations occurs in two distinct 
phases. First, the cloud as a whole forms cores of size 
~ 20,000 AU (0.1 pc). This size scale is half the Jeans 



12 



length at density p. These large scale cores each have a 
major filament within them that is above the critical line 
density for stability (Larson 1985 Inutsuka & Miyama 



1992, 1997), 



Ac 



2kT 
fimnG' 



(23) 



At our temperature, T 
is A cri t = 1-0 x 10 16 g cm 
filaments in our cores range 



10 K, the critical line density 
1 . The line densities of the 



from 1.7 x 10 16 to 4.0 x 



10 g c m , similar to fi lament line densities seen in 



Serpens ( Andre et al.||2010| ). The gene ral morpholo gy is 



norpi 

similar to the h ub-tilament structure in |Myers| (|2Q09[) and 
yers (2011a) with the hub at the center of eac 
The large scale fragmentation is not affecte d by radia- 



tion. In contrast to high-mass star formation flKrum holz 
et al. 2QQ7a| | Cunningham et aL 2011} Krumnolz et al. 



20lTf , tE ere is simply not enough protostellar luminos 



ity to affect the large scales except possibly at late times 
in the R simulation. The winds do travel through the 
entire simulation domain and could theoretically affect 
the fragmentation, but this does not happen in practice 
since the winds do not couple well to the cores. 

The second stage of fragmentation occurs as the fila- 
ment in each core contracts under self-gravity. The fil- 
aments can then fragment into many stars around the 
central stellar system. Unlike the large-scale core frag- 
mentation, this small scale fragmenta tion can be signif- 
icantly suppressed by radiation ( |Krumholz et al.||2007a| 
Offner et al.|2009b| ). This is best demonstrated in Figure 
7| Simulation K stagnates at 30 stars, while B and BW 
finish with 80 stars. At the end of R, the total protostel- 
lar luminosity is 2500 Lq and this is currently heating 
30 Mq of gas. Winds by themselves do not affect small 
scale fragmentation. The total number of fragments, and 
therefore the total number of stars, are the same in the 
barotropic simulations with and without winds, as shown 
in Figure [7| It should be noted that our fragments occur 
in the core and not in accretion disks. This agrees with 
highe r resolution core evolution simulations ( | Offner et al. 
|2010[ ). 

When winds and radiation are combined, the winds 
have a significant indirect effect. The winds lower the 
mass and accretion rate of the protostars, which low- 
ers the luminosity as seen in Figure [9] This means that 
the fragmentation suppression seen in comparing R to B 
should be reduced when comparing RW to BW. To inves- 
tigate this in more detail, we have shown the mass func- 
tions of RW and BW, both at t = 0.83%, in Figure [16. 
For purposes of comparison, we have excluded stars with 
mass M < O.O5M , since these stars are usually short 
lived and possibly subject to details of numerical sink 
particle formation or merging. The two mass functions 
are nearly the same. A two sample Kolmogorov-Smirnov 
test gives a 50% chance that both samples are drawn 
from the same underlying function. The primary differ- 
ence between the two mass functions lies in the heaviest 
star in the simulation. The core that forms the heaviest 
star fragments early in BW, turning a star that is 2.8 
M in RW into a 2.0 M star in BW with two extra 
M ~ O.4M stars. The reduced fragmentation in RW 
for that star might be expected from radiative feedback, 
since that is the most luminous star in the simulation. 

To understand why radiation does not significantly al- 




FlG. 16. — Mass function for all stars with mass > O.O5M0 for 
simulations RW (top) and BW (bottom). The Chabrier IMF is 
plotted to guide the eye, but the simulations are still accreting and 
are not expected to match the IMF. The radiative and barotropic 
mass functions are similar, showing that radiation does not signif- 
icantly affect the IMF in regions of low-mass star formation when 
winds are present. The main difference comes from the largest star 
in RW fragmenting more in BW, lowering the largest mass and 
adding more stars around 0.5Mq. 

ter star formation for regions of low-mass star formation 
such as we are consid ering, we ex pand on a line of rea- 
soning developed by Bate (2009b). He introduced an 
effective Jeans length, which we write as A e ff, defined by 
the condition 

GM 

A e ff = — , 24 

so that A e ff is the radius at which the escape velocity 
equals the isothermal sound speed, c s . Since the mass 
inside A e ff is M = (47r/3)pA 3 ff , one finds that Bates' ef- 
fective Jeans length is smaller than the standard Jeans 
length, Aj, which is given in equation (|3|, by a factor 
7r(4/3) 1//2 = 3.63. On the other hand, A e ff is comparable 
to the radius of the fragments we find in our simulations, 
which have a diameter of Aj/2, corresponding to a radius 
of Aj/4. Bates' effective Jeans mass, which is the mass 
within a radius A e ff, is smaller than the standard value 
in equation Q, by a factor 7r 3 / \/27 ~ 6.0. 

To estimate the dust temperature at a distance A e ff 
from a star of luminosity L, we follow [Bate (2009b) and 
ignore the frequency dependence of the dust absorption 
coefficient; A e ff is then given by the condition 



L = 4Tr\ 2 ef{ a SB T 4 



(25) 



where <jsb is the Stefan-Boltzmann constant. Since c s 
(/cT/zimn) 1 ^ 2 , we have 



1/5 



= 5.3 



io- 



-19 



g cm 



L 



1/5 



K. 



V 3/cctsb , 

(26) 

(Bate|l2009bj) adopted a fiducial luminosity of 150 L©, 
which gives T ~ 15 K for p c± 10 -19 g cm -3 . This is 
larger than the background temperature of 10 K that 
we have assumed; the Jeans mass is therefore increased 
and fragmentation is suppressed. However, the observed 
median luminosity in local star-forming regions is far 
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smaller: Enoch et al. (2009) find a median luminosity of 
only 1.5 Lq, which gives a temperature of only 5.8 K at 
A e ff. Since this is less than the background temperature, 
we conclude that protostars typically do not suppress 
fragmentation at densities of order 10 -19 g cm -3 , which 
are characteristic of low-mass star-forming regions. At 
higher densities, the accretion luminosity can raise the 
temperature enough to begin to suppress fragmentation 
in clouds with T ~ 10 K. 

4.5. Fragmentation in Rho Ophiuchus 

Our cloud fragments at about half the Jeans length 
(i.e. at ~0.1 pc), but then continues to fragment be- 
low this point. Frag mentation at the Jeans length is 
comm only observed ( [Blitz fc Williams|1997||Enoch et ah] 
2008). In instances where observers have the resolution 
and sensitivity to resolve fragmentation at scales below 
the Je ans length, however, even more fragmentation is 
found jMotte et al.||1998| jjohnstone et al.|[2QQQ| |Teix-| 
eira et al.||2UU7| |<Jnen fc Arce||2UlU| |Bontemps et al.| 



2010p . Fragmentation can be quantified as a function 
ol scale, r. Given a set of clump locations in a cloud, 
one can calculate a set of clump pair separations. Let 
the differential number of pairs separated by distance r 
be cL?Vp a i r = H(r)d\nr. The number of clump pair sepa- 
rations for randomly distributed clumps is H ran (r). The 
two point correla tion function, w can b e calculated from 
these quantities ( Johnstone et al.||2000 ), 



w(r) = 



H{r) 
H Tan (r) 



1. 



(27) 



The two-point correlation function has been measured 
for the central parsec of p Ophiuchus by Johnstone 
( |2000[ ). In this measurement, excess fragmenta- 

3 x 10 4 AU, similar to 



et al. 

tion (w > 0) is found below r 



the Jeans length of the cloud. There is a pow er law fit, 



Larson 



(1995) also mea- 
bund a power law 



w(r) oc r -0-75 , in this regime, 
sured clustering of stars in Taurus and : 
fit with a break at 8000 AU, but attributed the break to 
binary stars. The separation between stars has had time 
to evolve since the initial fra gmentation, so we narrow 



Johnstone et al. (2000). 



our focus to comparisons with 

To compare our simulations to the observed w, we first 
created optically thin column density maps of our sim- 
ulations and convolved them with a Gaussian with a 
FWHM of 1600 AU The resolution was ch osen to be 
similar to that from Johnstone et al. (|2000|). We u sed 
the Clumpfind algorithm from Williams et al.| ( |1994|) on 
the convolved column density map from simulation RW 
to obtain a list of clumps and their positions. To in- 
vestigate the possibility of time evolution of w, this is 
performed at an early time in the simulation and then 
again at a late time, t ~ 0.4% and t ~ 0.75%, respec- 
tively. The results are shown in Figure 17 As expected, 
the correlation function drops off above r ~ 4 x 10 4 AU, 
about 2/3 the Jeans length at p for our simulation. More 
remarkably, the correlation function in our simulation 
matches that measured in p Ophiuchus quite well at all 
scales below this drop off. The early and late time sim- 
ulations also generally agree with each other, suggesting 
there is not much time evolution in w. There is a discrep- 
ancy between the two times at r ~ 5 x 10 3 AU. At these 
small scales, fragmentation does increase in time, as high 
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Fig. 17. — Two point correlation function, w, for clumps in sim- 
ulation RW convolved with a 1600 AU Gaussian beam. The corre- 
lation function is calculated at an early time, t ~ 0.4% and a late 
time, t ~ 0.75£ff . The correlation is similar at early and late times, 
except for the smallest scales, where fragmentation increases over 
time. For comparison, the fit to r < 3 x 10 4 AU measurements 
from p Ophiuchus is also included. 

density regions have more time to form and fragment. 

It should be noted that at t ~ 0.4%, our stars have 
not provided very much feedback, and the simulation can 
be considered solely gravito- hydro dynamic. Simulations 
with just hydrodynamics and gravity are scale-free. The 
exact match with p Ophiuchus is partially due to our 
choice of cloud parameters to mimic p Ophiuchus. The 
scale- free results are the break in w(r) at the Jeans length 
the w(r) oc r -0,75 functional form. 

4.6. Observed Core Mass Functions 

There is a wealth of observat ions cataloguing masses 
of cores in star- forming regions ( [ Motte et al.||1998| Testi 
I fc Sargent|1998)|Johnstone et al.|2UUU[|2UUl|IMotte et al. 
20011 JBeuther fc S chilke 2004 ; IStanke et al.|2006| |Alves 
[et al.|2007||Enoch et al.|2008[|Sadavoy et al |2010p . Most 
ol these observations are unable to resolve the small scale 
fragmentation seen in p Ophiuchus, but still provide valu- 
able information. To recreate the observations, we pro- 
duced optically thin column density maps of our simu- 
lations in all three directions and convolved them with 
Gaussian beams chosen to match the observations. We 
then applied Clumpfind to the convolved column density, 
similar to our comparison to p Ophiuchus. 

The observations have a wide range of beam sizes due 
to the range in distances to star-forming regions, so we 
also used a range of beam sizes for comparison. We find 
that the CMF derived from our simulated observations is 
highly sensitive to the beam size used. As the beam size 
increases, the smallest cores are no longer detectable and 
drop out of the CMF. In addition, tight clusters of cores 
become unresolved and look like new, much larger cores. 
Both effects increase the median clump mass. T he effect 
of overlapping cores has been explored further in |Kainu~ 



lainen et al.|Q2009bJ ). This sensitivity of the CMF to res- 
olution is seen also in the observations, and shown in Fig- 



ure 



18 To gather the median m ass of a range of CMFs, 
we used the data tabulated in |Reid fc Wilson| p006). 
All cores from |Reid fc Wi lson ( 2006]) are detected using 
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Fig. 18. — Median mass of the CMF found in a cloud as a func- 
tion of resolution of the observation. The CMFs from synthetic 
observations of RW are green tria ngles. For comparison, CMFs 
tabulated in Reid &; Wilson ( 2006|> are included as well as 3 CMFs 
from [E noch et al. (2UU8J. The two lowest points from [Reid &:| 
|Wilson| ( |2UUtj| ) are p Uphiuchus at different wavelengths. 



Clumpfind. To supplement that data, we also used the 
three clouds (Serp e ns, Perseus and p Ophiuchus) from 



Enoch et al. 



(2008) 



Enoch et al 
on Clumpmid, but use a methoc 

suits. Clumpfind has many limitations (Pineda et al 



( |2008D do not rely solely 
that returns similar re- 



2009 Goodman et al. 20091), but is still useful lor the 



purposes of comparison. When compared to these ob- 
servational data, our simulated CMFs match quite well. 
The simulated CMF with 3200 AU and 6400 AU beams 
are interesting in particular, as this is the usual range 
of telescope beams for nearby clouds. The extent of our 
simulation domain (130,000 AU) prevents useful compar- 
isons to more distant observations. In addition, our base 
grid resolution, 512 AU, causes discrepancies with obser- 
vations of cores at small beam sizes. Bound cores will 
trigger adaptive refinement and can go to smaller scales, 
but smaller, ephemeral, cores that would show up in ob- 
servations are not always captured in the simulations. 
This means the simulated median masses at 800 AU and 
1600 AU are too high compared to observations which 
include all cores, bound or not. This is most apparent 
when comparing to the observations of p Ophiuchus with 
a beam size of 1600 AU. The net effect is a flattening in 
the clump mass versus beam size relation at small beam 
sizes for the simulated clump observations. 
Figure [18] does not include the CMF from the Pipe 
"Wes et al.|2007| ). 



Nebula 



The resolution is comparable 
to the best observations of p Ophiuchus, but the median 
mass is near 1 M@. This makes it stand apart from the 
other observations. This region has a low density and 
is not actively forming sta rs. In addition, th e observed 
cores are largely unbound ( Lada et "al7||2008 ). The Pipe 
Nebula should perhaps be considered a pre-star-forming 
region, unlike Serpens or p Ophiuchus. These pre-star- 
forming regions have much higher Jeans lengths and 
masses and have funda mentally different column density 
distribution functions ( Kainulainen et aL]|2009a ). 



Given that we have identified clumps at the beginning 
of star formation and have stars at the end of the simula- 
tion, a natural question is how the initial clumps relate to 
the final stars. We will first consider the cores identified 
by clumpfind discussed in the previous section. These 
cores represent what an observer might find in a sim- 
ilar cloud. For this comparison, we need to use simu- 
lations that have run to completion, so cannot use the 
radiative simulations. As a first pass at the correlation 
between the CMF and IMF, we will start with simula- 
tion B, where winds are absent and one might expect the 
CMF and IMF to overlap. The initial CMF and the final 
IMF for B are shown in Figure 19 For proper compar- 
ison, the mass functions are not normalized, and total 
counts at each mass are shown. To avoid triple counting 
cores in the CMF, the synthetic observations are taken 
in only one line-of-sight direction instead of all three. 
Otherwise these CMFs are the same as those in previ- 
ous sections. When focusing on the CMF found with 




4.7. CMF to IMF relation 
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Fig. 19. — Mass functions for both the initial cores found using 
clumpfind and for the final masses in stars. The y axis represents 
total counts and is not normalized. Top panel: cores found using 
a 1600 AU beam size. Bottom panel: cores found using an 800 AU 
beam size. 

the 1600 AU beam size, the initial CMF and final IMF 
are well correlated. The typical mass and total number 
of objects match well between cores and stars. Unfortu- 
nately, this comparison only holds for the 1600 AU beam 
size. The CMF is highly sensitive to beam size, as shown 
in Figure 18 When a slightly smaller, 800 AU, beam 
size is usedTthe clumps are too small and too numerous 
to all correlate with stars. It should be noted that this 
analysis is performed using a single clump identification 
method (Clumpfind) on a crowded field. It is possible 
that different methods of clump identification or clouds 
with well-separated clumps do not have as much sensi- 
tivity to telescope resolution in the resulting CMF. Even 
without resolution effects, we are considering all cores, 
bound and unbound. Lo w-mass cores are less likely to 
be bound or to form stars ( |Myers|2009| jPadoan fc Nord- 
lund 2011), which will skew the CMF to lower masses. 
CMFs that only include bound cores should be more cor- 
related with the IMF. 
We do not need to limit ourselves to synthetic obser- 
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vations of our simulations and can identify cores from 
the full 6 dimensional phase space. This is done us- 
i ng the ' find clump s' routine in the yt analysis toolkit 



( [Turk et al . 
more detail 



2QlH). 



Smith et al. 



contours to return a hierarc. 



The algorithm is described m 
" ( |20Q9D . It uses density 
ly ol clumps, where each 
clump can contain smaller child clumps. Our simula- 
tions have hundreds of local maxima in density large 
enough to be considered clumps. At t = 0.3tff, only 
40 of these are bound, defined as | potential energy | > 
(thermal energy + kinetic energy). We will only con- 
sider the bound clumps in this analysis. 

The bound clumps cover a large range of sizes with 
the largest clump filling almost the entire box (181 M©). 
This clump contains five child clumps. Each of these 
clumps contain many generations of children and grand- 
children. When comparing to the IMF, one should ig- 
nore clumps with multiple child clumps, and count the 
children instead. Throwing out clumps with masses less 
than O.O5M , there are only 5 bound childless clumps 
at t = 0.3tff, for a total mass of 2.5 M . At t = 0.4tff , 
some of the parent clumps split into more children, re- 
sulting in 16 bound childless clumps, though the total 
mass is largely unchanged, at 2.6 M . These clumps 
are notably smaller than the primary turbulent cores de- 
scribed in §3.1 and will be called 'sub-cores'. These sub- 
cores lead to the burst of star formation at t = OMf? 
and their mass function is shown in Figure 20 The me- 




-1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 

log(Mass/M Q > 

Fig. 20. — Histogram of masses of bound childless sub-cores at 
t = OAtff for simulation RW. 



dian mass of these sub-cores is 0.13 M . Even if each 
sub-core forms exactly 1 star, they do not explain the 
stellar mass function of the simulation, which eventually 
has a median mass of 0.3 M when winds and radia- 
tion are not included. This discrepancy can be explained 
by the 20,000 AU turbulent cores. These objects are 
bound, but they cannot be detected with density con- 
tours due to their supersonic turbulence. Each core has 
multiple pockets of high density gas. When using den- 
sity contours, one sees many smaller unrelated, unbound 
clumps instead of a few larger bound clumps that cor- 
respond to the physical cores. The cores are visible by 
eye and should be indentifiable with a more advanced 



density search. We identify them by stellar clustering. 
The central stellar systems (usually binary systems) in 
these cores are all much more massive than the sub-cores. 
If these central stellar systems are ignored, the median 
mass of stars moves from 0.3 M to 0.16 M , much 
closer to the median sub-core mass. The median mass 
of non-central stars for the case with winds is 0.07 M©. 
The turbulent core and central star properties are sum- 
marized in table lU 

The cores are bound at late times even when only grav- 
ity due to the stars is considered, with the kinetic energy 
in stars approximately half the potential energy of stars 
in each core. In addition, the core-to-core velocity dis- 
persion is typically 0.4-0.5 km s _1 . This is notably lower 
than the cloud velocity dispersion, which starts at 1.2 
km s _1 . There is a strong anti-correlation between veloc- 
ity and density, as the densest gas occurs at stagnation 
points in a turbulent flow. This means that the core- 
to-core velocity dispersio n will naturally be much lower 
than that of the gas flOffner et aD|2009a| . 

Once the cores have formed, each core is carved out by 
the outflows of its own protostellar system. This yields 
the core to star efficiency factor, 0.2 < e core < 1.0. The 
amount of mass lost from a spheroidal core can be cal- 
culated from the to tal momentum output and opening 
angle of the winds (Matzner & McKee 2000), but the 
cores in our simulations are more complicated. 

The best way to calculate e cor e is to compare the mass 
of stars in simulations with and without winds. In the 
simulations with the barotropic equation of state (B and 
BW), the total, the mean, and the median mass of stars 
are all approximately 3:1 comparing the non-wind simu- 
lation to the wind simulation at any point in time. This 
means e core ~ 1/3. The 2/3 of the core mass that is 
lost in the presence of winds is mostly core gas that is 
entrained in the winds. This was not precisely quan- 
tified, and it is possible other mechanisms account for 
some mass loss, such as unbinding outer core gas when 
the mass from the interior core is lost. In the radia- 
tive simulations (R and RW), the total mass of stars is 
also 3:1 comparing the two simulations. The mean and 
median masses are 3:1 for simulation RW, but closer to 
10:1 in simulation R due to unphysically strong radiative 
suppression of fragmentation. 

In the event that cores are accreting mass, e core is a 
function of time, and may be larger than 1 if the instanta- 
neous core mass is less than the final stellar system mass 
(Padoan fc Nordlund||2011| ). Given that the sub-cores in 



our simulations are massive enough at early times to cre- 
ate their final stellar systems, the sub-cores probably do 
not accrete much mass. It is extremely difficult to track 
these cores over time, however, so we cannot quantify the 
amount of sub-core accretion. 

4.8. Turbulent Fragmentation and Competitive 
Accretion 

It is useful to place the stellar accretion in our simu- 
lations in the context of existing star formation models. 
Two current popular models are turbulent fragmentation 
and competit ive accretion. In the turbulent fragmenta- 
tion model, flPadoan||1995l |Padoan fc Nordlund| [2 002 
McKee fc Tan||Mfl [Padoan et al.||M7| |HennebeiTe"T 



Chab rier||2QQ8[ |2009p , supersonic turbulence in molecu- 



lar clouds creates many cores. Each bound core then 
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collapses into a single stellar system (if the structure of 
the core is dominated by internal turbulence, the result- 
ing model for the formatio n of the star was ter med the 



"turbulent core model" by |McKee fc Tan||2QQ3| ). In this 
scenario, the mass from each star is accreted almost en 
tirely f rom its natal co re. In the competitive accretion 
model 



M 

mel. 



Zinnecker 1982; 



fc Bonriell||2005j IBonnell et al. 



begin with a mass 



Bonnell et al 
2006), 



1997, 2001 Bate 



the bound cores 



The molecular cloud 



0.1 Mq. 

undergoes a global collapse and all stars accrete from 
the entire cloud. Protostars exhaust their cores at low- 
masses and then grow by Bondi-Hoyle accretion. The 
protostars compete with each other for mass from the 
host cloud and the dynamics of this competition lead to 
the IMF. Roughly speaking, the virial parameter of the 
cloud decides which model is applicable ([Krumholz et al. l 



2005 Bonnell & Bate 2006 Offner et 



with sufficiently sub-virial turbulence, g. 



1012008). 
:e, global < 



In clouds 
collapse is 

possible and competitive accretion prevails. In virialized 
clouds, core accretion dominates. Most simulations of 
star formation start with virialized clouds, but turbu- 
lence quickly dissipates and simulations that do not re- 
generate turbulence become sub-virial and demonstrate 
competitive accretion. Turbulence can be regenerated 
by protostellar outflows, HII regions, or a cascade from 
larger scales; in simulations, the cascade can be regener- 
ated by large scale driving. 

Our simulations largely agree with turbulent fragmen- 
tation models, while introducing a hierarchical aspect 
of sub-cores within turbulent cores. We form turbulent 
cores on the cloud Jeans length and each of those cores 
forms a central binary or single star with mass roughly 
equal to the core mass. Even the smaller stars are formed 
in their own sub-cores and do not accrete from the cloud 
at large. Our turbulent cores do accrete from the larger 
cloud (increasing their initial mass by ~ 75% over the 
course of the simulations), which was not originally part 
of the core accretion theory, but recent turbulence simu- 
lations suj | rges^ their host 
cloud ( jFalceta-Goncalves fc Lazarian||201l| ) . 

The accretion from the larger cloud onto cores in our 
simulation is possibly caused by the fact that we have 
not included turbulent driving from large scales. Turbu- 
lence is regenerated to some extent by protostellar out- 
flows, but this is relatively ineffective in denser gas. Our 
simulations are then similar to the undriven simulations 
which also show accretion onto 



m 



cores. 



Offn er et al. 
The simu" 



(|2008j) 

at ions in 



Offner et al. (2008) with ex- 



ternal driving produce cores that do not accrete much 
from the cloud. If magnetic fields were included in our 
simulati ons, the outflows would couple to much more of 
the gas flWang et al.|2010 ), which would reduce accretion 
onto the cores. 

4.9. The Role of Stellar Mergers 

The stellar mass functions in our simulations are in- 
fluenced by the details of our sink particle merger pro- 
cess. Mergers are necessary because all codes will intro- 
duce numerical fragmentation once they c an no longer re- 



solve the Jeans length on the finest scale (Truelove et al. 
1997). More stringent sink particle conditio ns can re^ 



duce" the n umber of unwanted sink particles (Federrath 
et al.|2010| ), but numerical fragments are unavoidable. 11 



fragments will steal mass from physical fragments and 
masquerade as real stars, artificially lowering the IMF. 
On the other hand, allowing sink particles to merge has 
a similar effect on the IMF as suppressing gas fragmen- 
tation. If all sink particles that pass near each other are 
merged together, the particles will consolidate over time. 
Eventually, all IMFs look similar to the top heavy IMF 
from R, where radiation suppressed most fragmentation. 
Our decision to only merge protostars with masses less 
than O.O5M , the mass of the first core of a 1 M star, 
is a compromise between the two extremes of no mergers 
and all mergers. 

This suppression of mergers is seen comparing the 
IMFs of RW and BW. The radiative simulation should 
suppress some fragmentation while the barotropic sim- 
ulation fragments all the way down to the Jeans mass 
at the self-opacity limit, m ~ O.OO4M . Nevertheless, 
they have the same IMF. The barotropic simulation does 
fragment much more than the radiative simulation near 
the resolution limit. This is seen in the total number of 
sink particles created, where the barotropic simulation 
creates 7 times more particles than the radiative simula- 
tion. These particles are nearly all very small in mass and 
immediately merged. The net effect of the extra merg- 
ers is to suppress fragmentation by combining fragments 
below the merger radius (128 AU). The two simulations 
are nearly identical above this scale and therefore pro- 
duce the same IMF. 

Our choice of merger mass (m merge = 0.05 M ) is 
based on calculations of first collapse of a solar mass 



star ( Masunaga et al.|1998 Masunaga fc Inutsuka| 2000 



but the correct mass is not certain. In addition, SPH 
simulations with resolutions of a few AU (compared to 
our 32 AU) move their numerical fragmentation to much 
small er scales and do not show stella r mergers o f this 
type ( IStamatellos fc Whitworth||2009| |Bate|j2011| ). To 
investigate the effect of our mass choice, we repeated 
RW with a merger mass of 0.01 M out to t = 0.55 tff. 
At this point in the simulations, the total mass in stars is 
2 Mq. In the simulation with m merge = 0.05 M , there 
are 17 total sink particles; whereas in the simulation with 
Emerge = 0.01 M , there are 30 total sink particles, even 
though the total mass in particles is the same. Given that 



the final IMF for 



Emerge 



0.05 Mq is at slightly lower 



masses than the observed IMF, nearly doubling the num- 
ber of sink particles with m merge = 0.01 M would move 
the IMF to masses less than half the observed value. 



sink particles are not allowed to merge, these numerical 



4.10. Scaling 

For most isothermal simulations, there is a free param- 
eter which allows rescaling of the mass, M, density, p, or 
temperature T, while maintaining all of the same dimen- 
sionless parameters a, .M, and Mj/M. When the merger 
mass is included, m merge also scales with M. Now if we 
scale the simulation to a higher mass, we are also increas- 
ing m mer ge- There is some uncertainty in m mer ge, but it 
would be difficult to justify increasing it much more than 
our current level. Even increasing m merge by a factor of 
2 would bring it to uncomfortably close to the character- 
istic IMF mass of 0.2 M . This means the mass scales of 
our simulation are relatively stationary. When protostel- 
lar winds are included, they introduce a new fixed dimen- 
sionless number, the Mach number of the winds .M w ind- 
Since the speed of the winds is proportional to the escape 
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velocity from the stellar surface (i.e., oc M 1 / 2 [McKee fc 
Ostriker||2007|), A^ w i n d sets the quantity M/T. In prac- 



tice, this is not a very tight co nstraint because the wind 
speed itself is quite uncertain ( |Downes fc Cabrit]|2007 ). 
When radiation is important, the luminosity ol each star 
is set by complicated stellar models that depend on M 
as well as the accretion history. The time scale, which 
goes into the accretion rate, is set by % and therefore p. 
In addition, the resulting radiation-hydrodynamics de- 
pends on the temperature. This uniquely sets all scales 
in the simulations. Even when radiation is not dynam- 
ically important, we do match the observed protostel- 
lar luminosities and cannot change our masses without 
jeopardizing the agreement. Using the approximation 
L oc our cloud mass is constrained to 

160 M Q < M < 200 M before it no longer falls in the 
error bars luminosities. Cloud masses below the lower 
limit do not match the observed median luminosity and 
masses above the upper limit do not match the mean. 
Our general match to the IMF, which sets M, provides 
an additional, looser constraint. 

5. CONCLUSIONS 

We report the results of several simulations of the for- 
mation of a low-mass star-forming cluster, comparable 
to the central parsec of p Ophiuchus. Our simulations 
achieve 32 AU resolution using adaptive mesh refine- 
ment. We also include radiation- hydrodynamics and pro- 
tostellar feedback. The protostellar feedback represents 
both protostellar radiation and bipolar outflows. To iso- 
late the individual effects of radiation and outflow feed- 
back, we perform a suite of 4 simulations: a base simu- 
lation with no feedback (i.e., with a barotropic equation 
of state and no outflows), a radiative simulation with no 
outflows, an outflow simulation with no radiation and a 
simulation with both outflows and radiation. This is the 
first simulation of a low-mass star-forming cluster with 
both radiation and protostellar outflows. 

When outflows are included in the simulation, they are 
able to replenish the kinetic energy lost from decaying 
turbulence. This new outflow-driven turbulence is fun- 
damentally different than isotropi c turbulence driven on 
large scales, however, as seen in (Nakamura & Li||2007). 



This can be measured in the solenoidal to compressional 
energy ratio, which climbs from 2 to 7 over the course 
of our simulations. Hydrodynamic outflow-driven turbu- 
lence does not couple well to the high density gas and 
cannot prevent cores from collapsing. 

Both simulations with outflows reproduce the expected 
mass functions. The radiative simulation does not finish 
accreting, but it matches t he turbulent core proto stellar 
mass function (PMF) from |McKee fc Offner | ( [gOlQj ) . The 
barotropic si mulation has mo stly finished accretion and 
matches the |Chabrier| ( 2005 ) IMF. Simulations without 
winds produce mass functions that are too massive by 
factors of 3 and 10 for the barotropic and radiative sim- 
ulations, respectively. 

When we compare final stellar masses with and with- 
out outflows, we find the outflows remove 2/3 of the mass 
that would go into stars. This creates a core effic iency pa- 
rameter e CO re — 1 /3, similar to predictions from Matzner 
( |2Q00| . 



& McKee 



The importance of outflows in our 
IMF calls into question simulations that produce the ob- 
served IMF without outflows (e.g. Bate 2009b; Price & 



Bate 2009). The outflows do not significantly affect the 
overall cloud dynamics, as they have small opening an- 
gles and do not couple well to the dense gas in the cores. 
It is likely magnet ic fields would change that conclusion 
flWang et al.]|201(^ . 

Outflows also significantly reduce protostellar lumi- 
nosities. They reduce the accretion rate onto a protostar 
by a factor of 3 by disrupting its parent core, and there- 
fore reduce the final mass of the protostar by a factor 
of 3. The typical radius of a protostar does not change 
much in this mass regime, so accretion luminosity is re- 
duced by an order of magnitude. This luminosity lost 
this way is much larger than the mechanical luminosity 
of the outflows themselves. 

The reduced luminosities in the simulation with radia- 
tion and outflows allows it to match the observed proto- 
stellar luminosities of nearby star- forming regions. This 
includes the mean and median luminosities as well as 
the standard deviation of the log luminosity. The proto- 
stellar luminosity function also depends on the tim e and 
mass evolution of the protostellar accretion rate( Offner 



fc McKee 2011). We find the accretion rate scales with 
mass as predicted from core accretion model s that in- 



clude both thermal and nonthe rmal motions (Myers 
Fuller|1992[ |McKee fc Tan|2003[ ) as opposed to competi- 



tive accretion or isothermal sphere models. We also find 
the accretion rate mu st be tapered and is cons istent with 
the linear tapering in |McKee fc Offner] ( |2 01 0\ . Since our 
models do not resolve disk physics, they do not include 
episodic accretion. The fact that we nonetheless are able 
to reproduce the observed protostellar luminosities in re- 
gions of low-mass star formation suggests that episodic 
accretion is not the dominant factor in resolving the lumi- 
nosity problem for low -mass protostars; this is c onsistent 
with the dis cussion of |Offner ~fc McKee] (2011 ), but not 
with that of |Dunham et al.| ( |2(J10| ). 

The simulation with radiation and without outflows 



confirms the finding of Offner et al. (2009b) and 
|Krumholz et al. (2011) that protostars can heat their 
host cloud enough to suppress fragmentation. When out- 
flows are included, however, the total luminosity of stars 
drops by an order of magnitude, and radiation is far less 
effective at suppressing fragmentation. The simulation 
with radiation and winds has over twice as many stars as 
the simulation with radiation but without winds. When 
fragmentation is additionally suppressed by merging pre- 
second-collapse stars at 128 AU, radiation has almost no 
effect on the resulting mass function. Thus, radiation 
may be necessary to suppress fragmentation below ~128 
AU for the conditions simulated here, but it does not 
significantly effect the gas dynamics above those scales. 

We are able to recreate the clustering properties of the 
cores found in p Ophiuchus, measured by the two point 
correlation function, w(r) of observed clumps. Our sim- 
ulated observations recreate the overall normalization, 
the w(r) oc r -3 / 4 slope below the Jeans length, and the 
break in the power law at the Jeans length observed by 
| Johnstone et al.| ( |2000[ ). This implies our simulation of 
fragmentation is accurate down to at least 2,000 AU, set 
by the resolution of the observations. 

To investigate the conversion of the observed core mass 
function (CMF) to the stellar IMF, we create simulated 
dust maps and find cores using Clumpfind. The mean 
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core mass in observations systematically increases with 
telescope resolution, and we are able to recreate those 
core masses over an order of magnitude of resolutions. At 
resolutions typical of observations of nearby star-forming 
regions, 1600 AU, the early time CMF and the final 
IMF overlap when outflows are not included. Outflows 
shift the IMF to lower masses while retaining the orig- 
inal shape. When we simulate observations with better 
resolutions, however, the CMF changes and no longer 
matches the IMF. This suggests that current observa- 
tions suggesting a CMF to IMF correspondence could 
change at higher resolutions. It should be noted that 
this analysis was performed with Clumpfind on a simu- 
lated cloud with closely spaced cores. It is possible that 
different methods of clump identification or a cloud with 
more separation between cores would produce a CMF 
that is not as sensitive to telescope resolution. 

In sum, we have used ORION, an adaptive mesh re- 
finement (AMR) gravito-radiation-hydrodynamics code, 
to simulate low-mass star formation in a turbulent molec- 
ular cloud in the presence of protostellar feedback. Our 
results for the most realistic simulation, which includes 
both radiation and feedback, are consistent with obser- 
vations of the protostellar luminosity function, the core 
mass function, the two-point correlation function of cores 
in p Ophiuchus, and the protostellar and initial stellar 
mass functions. We find that protostellar radiation does 
not affect fragmentation below 128 AU and that pro- 
tostellar outflows do not significantly affect large-scale 
cloud dynamics, serving only to reduce the mass in cores 
for a core to star efficiency, e cor e — 1/3. Lastly, we 



find the accretion histories of our stars match core ac- 
cretion models that include both thermal and turbulent 
motions and appear to be inconsistent with competitive 
accretion or isothermal sphere accretion. The inclusion 
of magnetic fields would increase the coupling between 
outflows and dense gas, possibly changing the large-scale 
cloud dynamics. With an increase in resolution, it would 
be possible to explicitly model protostellar mergers prior 
to second collapse, as well as investigate the small-scale 
regime where protostellar radiation suppresses fragmen- 
tation. 
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TABLE 3 

Turbulent Core Properties for Barotropic Simulations. 





Without Winds 


With Winds 


Mg as Initial 


M gas Final 


M*, total 


M* ? central 


M g as Final 


M* t total 


-M* ? central 


3.4 


3.6 


5.6 


4.7 


2.3 


1.9 


1.3 


7.0 


3.1 a 


16.8 a 


9.4 


2.0 a 


9.6 a 


5.9 


3.8 


2.5 


5.0 


3.4 


1.2 


1.5 


1.4 


5.0 


2.7 


5.1 


3.9 


1.6 


2.9 


2.2 


2.9 


_a 


_a 


1.0 


_a 


_a 


0.04 


3.3 


_a 


_a 


2.7 


_a 


_a 


0.3 



a The last two cores eventually merge with the largest core, making 
it impossible to measure the final gas and total stellar properties. 
The final properties of the largest core are necessarily a sum over 
the last two cores in addition to the largest core. 



